﻿/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms which obtain the physical parameters of the Sun
Created: PJN / 29-12-2003
History: PJN / 16-06-2004 1) Fixed a typo in the calculation of SunLongDash in Calculate.
                          Thanks to Brian Orme for spotting this problem.

Copyright (c) 2003 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/


using System;

namespace AAPlus
{
    public class PhysicalSunDetails
    {
        public double P = 0.0;
        public double B0 = 0.0;
        public double L0 = 0.0;

        public PhysicalSunDetails()
        {
        }
    }

    public static class PhysicalSun
    {

        public static PhysicalSunDetails Calculate(double JD)
        {
            double theta = CoordinateTransformation.MapTo0To360Range((JD - 2398220) * 360 / 25.38);
            double I = 7.25;
            double K = 73.6667 + 1.3958333 * (JD - 2396758) / 36525;

            //Calculate the apparent longitude of the sun (excluding the effect of nutation)
            double L = Earth.EclipticLongitude(JD);
            double R = Earth.RadiusVector(JD);
            double SunLong = L + 180 - CoordinateTransformation.DMSToDegrees(0, 0, 20.4898 / R);
            double SunLongDash = SunLong + CoordinateTransformation.DMSToDegrees(0, 0, Nutation.NutationInLongitude(JD));

            double epsilon = Nutation.TrueObliquityOfEcliptic(JD);

            //Convert to radians
            epsilon = CoordinateTransformation.DegreesToRadians(epsilon);
            SunLong = CoordinateTransformation.DegreesToRadians(SunLong);
            SunLongDash = CoordinateTransformation.DegreesToRadians(SunLongDash);
            K = CoordinateTransformation.DegreesToRadians(K);
            I = CoordinateTransformation.DegreesToRadians(I);
            theta = CoordinateTransformation.DegreesToRadians(theta);

            double x = Math.Atan(-Math.Cos(SunLong) * Math.Tan(epsilon));
            double y = Math.Atan(-Math.Cos(SunLong - K) * Math.Tan(I));

            PhysicalSunDetails details = new PhysicalSunDetails();

            details.P = CoordinateTransformation.RadiansToDegrees(x + y);
            details.B0 = CoordinateTransformation.RadiansToDegrees(Math.Asin(Math.Sin(SunLong - K) * Math.Sin(I)));

            double eta = Math.Atan(Math.Tan(SunLong - K) * Math.Cos(I));
            details.L0 = CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(eta - theta));

            return details;
        }

        public static double TimeOfStartOfRotation(long C)
        {
            double JED = 2398140.2270 + 27.2752316 * C;

            double M = CoordinateTransformation.MapTo0To360Range(281.96 + 26.882476 * C);
            M = CoordinateTransformation.DegreesToRadians(M);

            JED += (0.1454 * Math.Sin(M) - 0.0085 * Math.Sin(2 * M) - 0.0141 * Math.Cos(2 * M));

            return JED;
        }

    }
}
